Complications in Bronchoscopic Lung Volume Reduction

Insights form the National Inpatient Sample

Author

Hala Nas

Published

2024-04-28

R Setup

Code
knitr::opts_chunk$set(comment = NA) 

library(janitor)
library(broom)
library(haven)
library(naniar)
library(mice)
library(kableExtra)
library(tableone)
library(skimr)
library(Matching)
library(cobalt)
library(survival)
library(twang)
library(survey)
library(conflicted)
library(knitr)
library(tidyverse)

conflicts_prefer(dplyr::select(), dplyr::filter())

theme_set(theme_bw()) 

1 Background

Emphysema represents distortion in the lung parenchyma that occurs in patients with chronic obstructive pulmonary disease (COPD) and is pathologically defined as permanent enlargement of the alveoli along with destruction of the alveolar walls. Emphysema can lead to lung hyperinflation and air trapping, which contributes to refractory dyspnea. In carefully selected patients with emphysema, lung volume reduction could improve breathlessness, lung function, quality of life and exercise capacity by decreasing lung hyperinflation 1.

Two mechanisms exist for lung volume reduction, the first, lung volume reduction surgery (LVRS), entails wedge resection of the emphysematous tissue 2. Surgical morbidity can be significant, and certain comorbidities can preclude treatment in many patients. The second is bronchoscopic lung volume reduction (BLVR) with endobronchial valves (EBV). It is a procedure where one-way valves are placed into airway segments in the most diseased lung lobe, allowing air to escape the target area but not allowing it to go back in 3. This mechanism allows the eventual deflation of the target lung lobe and the offloading diaphragmatic stretch.

Two Large clinical trials 4,5 resulted in the FDA approval of endobronchial valves for bronchoscopic lung volume reduction. Overall, respiratory complications occur in 31–35% of patients 6. The most common adverse event associated with EBV placement is pneumothorax (approximately 25% of the cases) 7,8, particularly in the immediate post operative period. Other complications include respiratory failure, acute exacerbation of COPD (AECOPD), pneumonia and hemoptysis. The frequency of these complications is significantly lower than pneumothorax in the published clinical trials.

Since FDA approval in 2018, several centers in the United States have been offering bronchoscopic lung volume reduction for patient with severe COPD. Some centers have published individual reports of complication rates 9,10. To date, there has been no study evaluating complication rates post bronchoscopic lung volume reduction on a population level. While pneumothorax is an expected complication of BLVR, acute respiratory failure can lead to dramatic instability and if progressive, can result in endobronchial valve removal.

2 Overview of the Intervention

2.1 Bronchoscopic Lung Volume Reduction (BLVR)

2.2 Surgical Lung Volume Reductuion (LVRS)

3 Objectives

To evaluate the odds of acute respiratory failure (ARF) and pneumothorax (PTX) in adult patients with a principal diagnosis of emphysema who are undergoing elective bronchoscopic lung volume reduction (BLVR), in comparison to patients with emphysema who are undergoing surgical lung volume reduction (LVRS) in the immediate post operative period

4 Research Questions

Is bronchoscopic lung volume reduction associated with a higher odds of developing post-procedural acute respiratory failure in adult patients with emphysema compared to lung volume reduction surgery?

Is bronchoscopic lung volume reduction associated with a higher odds of developing post-procedural pneumothorax in adult patients with emphysema compared to lung volume reduction surgery?

5 The Data

Data were obtained from the national inpatient sample (NIS); a US administrative database published under the healthcare utilization project (HCUP) by the Agency for Healthcare Research and Quality 11 through a data use agreement.

The NIS is the largest publicly available all-payer inpatient care database in the United States, containing data on more than seven million hospital stays from nearly 1,000 U.S. hospitals including non-federal public hospitals and academic medical centers. The data include discharge-level records on demographics, diagnoses, procedures, length of hospital stay (LOS), and resource utilization. The NIS uses the International Classification of Diseases, Clinical Modification 10 (ICD-10) to code for diagnoses and procedures.

The data were de-identified, not requiring an IRB approval

6 The Participants

Adult patients with a primary diagnosis of emphysema admitted electively for lung volume reduction (surgical or bronchoscopic).

The inclusion and exclusion criteria are aimed at ensuring that the patient was solely admitted for lung volume reduction, and that surgery or endobronchial valve placement was not performed for persistent air leak.

Inclusion criteria:

  • Adults (18 years and older)
  • Elective admission
  • Principal diagnosis of emphysema (I10-DX1)
  • Principal procedure for the admission is lung volume reduction, including bronchoscopic lung volume reduction (I10-PR1)
  • Principal procedure on the day of admission (procedure on day 0)

Exclusion criteria

  • Non-elective admissions
  • Presence of persistent air leak diagnosis (ICD10 codes)

The sample size consists of 518 adult patients with 253 (49%) patients undergoing bronchoscopic lung volume reduction and 265 patients undergoing surgical lung volume reduction (51%).

The study is a retrospective cohort study, two groups of patients with the same baseline disease (emphysema), one group is receiving the intervention of choice (BLVR), and followed during admission for the development of the outcomes after the intervention of choice.

7 The Intervention

Bronchoscopic lung volume reduction (BLVR) is the intervention of interest. 253 (49%) of the participants out of 518 in my sample have received BLVR and will be compared to surgical lung volume reduction (LVRS), in 265 (51%) of the participants.

The intervention (BLVR) will be defined as the first procedure (I10_PR1) received by the participant on day 0 of the elective admission, using ICD-10 procedure codes: 0BH38GZ, 0BH88GZ, 0BH48GZ, 0BH58GZ, 0BH68GZ, 0BH78GZ, 0BH8GZ, 0BH98GZ.

The comparison (LVRS) will be defined as the first procedure (I10_PR1) received by the participant on day 0 of the elective admission, using ICD-10 procedure codes: 0BBC0ZZ, 0BBC4ZZ, 0BBD0ZZ, 0BBD4ZZ, 0BBF0ZZ, 0BBF4ZZ, 0BBG0ZZ, 0BBG4ZZ, 0BBH0ZZ, 0BBH4ZZ, 0BBJ0ZZ, 0BBJ4ZZ, 0BBK0ZZ, 0BBK4ZZ, 0BBL0ZZ, 0BBL4ZZ, 0BBM0ZZ, 0BBM4ZZ

8 The Covariates used in the propensity score

The propensity score model used covariates that are related to patient characteristics (age, sex, race, presence of chronic respiratory failure, insurance carrier and classification of the median income for the patient’s zip code) and hospital characteristics (hospital size, location, teaching status, hospital ownership and year of the procedure)

The covariates were measured and recorded on admission and prior to the procedure.

9 The Outcomes

9.1 Primary outcome: Acute respiratory failure (ARF)

Respiratory Failure

Acute respiratory failure is a new (or increased) need for oxygen supplementation in order to maintain adequate oxygenation.

It will be defined using ICD-10 diagnosis codes J9600, J9601, J9602, J9620, J9621, J9622, J95821, J95822 and is a binary outcome (Yes/No)

9.2 Secondary outcome: Pneumothorax (PTX)

Pneumothorax

Pneumothorax is the abnormal collection of air around the lung leading to lung collapse, and possibly causing respiratory failure or in severe forms, instability and death.

This will be defined using ICD-10 diagnosis codes J930, J9311, J9312, J9383, J939, J95811, J95812 and is a binary outcome (Yes/No)

10 Importing and cleaning the data

10.1 Importing the data

The final analytical data was prepared using SAS. This included importing the national inpatient sample (NIS) into a statistical program that can handle the large size and number of original observations (roughly 21 million observations). SAS was used to properly create the variables of interest based on ICD-10 codes. A final dataset with the intended variables was then created for the final analyses.

The data was imported in the format of sas7bdat using haven package. zap_labels and clean_names were used to remove column labels and standardize names respectively into R format.

Code
blvr <- read_sas("blvr.sas7bdat", NULL) |> zap_label() |> clean_names() 
Code
glimpse(blvr)
Rows: 540
Columns: 18
$ age              <dbl> 70, 69, 73, 57, 67, 55, 44, 55, 65, 43, 81, 69, 63, 7…
$ elective         <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ female           <dbl> 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0,…
$ pay1             <dbl> 1, 3, 1, 1, 3, 1, 3, 3, 1, 3, 1, 1, 3, 1, 3, 2, 2, 3,…
$ race             <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 6, 1, 2, 1,…
$ year             <dbl> 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018,…
$ zipinc_qrtl      <dbl> 2, 2, 4, 2, 2, 4, 2, 4, 4, 4, 3, 4, 4, 3, 3, 3, 1, 1,…
$ hosp_bedsize     <dbl> 2, 3, 3, 3, 2, 2, 3, 3, 3, 2, 1, 3, 3, 3, 1, 1, 3, 2,…
$ hosp_locteach    <dbl> 3, 3, 3, 3, 3, 3, 1, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3,…
$ hosp_region      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ h_contrl         <dbl> 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2,…
$ emphysema        <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ ebv              <dbl> 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,…
$ lvrs             <dbl> 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,…
$ crf              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ ptx              <dbl> 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1,…
$ arf              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,…
$ volume_reduction <chr> "LVRS", "BLVR", "LVRS", "BLVR", "LVRS", "LVRS", "LVRS…

10.2 Cleaning the names and creating factor variables

The code below achieved the following steps

  • Filtering to include adult patients only (Over the age of 18)
  • Renamingebv to vol_reduction, this was intended to be the numeric 0/1 version of the intervention
  • Mutating character variables as factors and appropriately coding them
  • Re-leveling certain factors appropriately
  • The tibble will continue to include vol_reduction which is the intervention of choice as a numeric variable (0/1) for further analysis
Code
blvr <- blvr |>
  filter(age > 17) |> rename(vol_reduction = ebv) |> 
  mutate(sex = fct_recode(factor(female), "Male" = "0", "Female" = "1"),
         insurance = fct_recode(factor(pay1), "Medicare" = "1", "Medicaid" = "2", "Private" = "3", "Private" = "4", "Private" = "6"),
         insurance = fct_relevel(insurance, "Medicare", "Medicaid", "Private"),
         race = fct_recode(factor(race), "White" = "1", "Black" = "2", "Other" = "3", "Other" = "4", "Other" = "5", "Other" = "6"),
         hospital_size = fct_recode(factor(hosp_bedsize), "Small" = "1", "Medium" = "2", "Large" = "3"),
         teaching = fct_recode(factor(hosp_locteach), "Rural" = "1", "Urban nonteaching" = "2", "Urban teaching" = "3"),
         region = fct_recode(factor(hosp_region), "Northeast" = "1", "Midwest" = "2", "South" = "3", "West" = "4"),
         ownership = fct_recode(factor(h_contrl), "Government" = "1", "Private, non-profit" = "2", "Private, for-profit" = "3"),
         income = fct_recode(factor(zipinc_qrtl), "0-25th percentile" = "1", "26th to 50th percentile" = "2", "51st to 75th percentile" = "3", "76th to 100th percentile" = "4"),
         crf = fct_recode(factor(crf), "Yes" = "1", "No" = "0"),
         ptx_f = fct_recode(factor(ptx), "Yes" = "1", "No" = "0"),
         arf_f = fct_recode(factor(arf), "Yes" = "1", "No" = "0"),
         vol_reduction_f = as_factor(volume_reduction),
         vol_reduction_f = fct_relevel(vol_reduction_f, "LVRS"), 
         year = as_factor(year)) 

10.3 Final analytical tibble

The final analytical tibble included the covariates used in building the propensity score, as well as both numeric/ factor versions of the intervention and outcomes

Code
blvr <- blvr |> select(age, sex, race, crf, income, insurance, hospital_size, teaching, region, ownership, year, vol_reduction, vol_reduction_f, arf, arf_f, ptx, ptx_f)

10.4 Missing variables and Imputation

Code
blvr |> miss_var_summary()
# A tibble: 17 × 3
   variable        n_miss pct_miss
   <chr>            <int>    <dbl>
 1 race                14     2.70
 2 income               9     1.74
 3 age                  0     0   
 4 sex                  0     0   
 5 crf                  0     0   
 6 insurance            0     0   
 7 hospital_size        0     0   
 8 teaching             0     0   
 9 region               0     0   
10 ownership            0     0   
11 year                 0     0   
12 vol_reduction        0     0   
13 vol_reduction_f      0     0   
14 arf                  0     0   
15 arf_f                0     0   
16 ptx                  0     0   
17 ptx_f                0     0   

There were fourteen missing observations in the race variable, and nine observations in the income variable. These constituted 4.4% of the total data, and assumed to be missing at random (MAR).

Code
blvr |> tabyl(income)
                   income   n    percent valid_percent
        0-25th percentile 139 0.26833977     0.2730845
  26th to 50th percentile 133 0.25675676     0.2612967
  51st to 75th percentile 129 0.24903475     0.2534381
 76th to 100th percentile 108 0.20849421     0.2121807
                     <NA>   9 0.01737452            NA
Code
blvr |> tabyl(race)
  race   n    percent valid_percent
 White 421 0.81274131    0.83531746
 Black  46 0.08880309    0.09126984
 Other  37 0.07142857    0.07341270
  <NA>  14 0.02702703            NA

Given the small sample and the desire to keep as many observations as possible, single imputation using the mice package. Fifteen imputed data sets were created and one dataset was later used as the single imputation tibble.

Code
set.seed(123)
blvr_mice <- mice(blvr, m = 15, printFlag = FALSE)
Code
blvr <- complete(blvr_mice, 10) |> tibble()

Checking the variables with the missing values after imputation showed preserved ratios in each level in comparison to pre-imputation

Code
blvr |> tabyl(income)
                   income   n   percent
        0-25th percentile 142 0.2741313
  26th to 50th percentile 136 0.2625483
  51st to 75th percentile 131 0.2528958
 76th to 100th percentile 109 0.2104247
Code
blvr |> tabyl(race)
  race   n    percent
 White 432 0.83397683
 Black  48 0.09266409
 Other  38 0.07335907

10.5 Adding an ID row

An ID row was added to identify observations if needed

Code
blvr <- tibble::rowid_to_column(blvr, "ID"); blvr <- blvr |> mutate (ID = as.character(ID))

10.6 Saving final tibbe

Code
write.csv(blvr, "blvr.csv")

11 Table (1) And Codebook

11.1 Final Codebook

Variable Type Description Role
age Numeric Age in years Co-variate
sex Factor (w/2 levels) Sex of the subject: Male or Female Covariate
race Factor (w/3 levels) White, Black, Other Covariate
crf Factor (w/2 levels) Presence of chronic respiratory failure: Yes or No Covariate
income Factor (w/4 levels) Quartile classification of the estimated median household income of residents in the patient’s ZIP Code: 0-25th percentile, 26-50th percentile, 51-75th percentile and 76-00th percentile Covariate
insurance Factor (w/3 levels) Medicare, Medicaid or Private Covariate
hospital_size Factor (w/3 levels) Small, Medium, Large Covariate
teaching Factor (w/3 levels) Teaching status of the hospital: Rural, Urban non-teaching and Urban teaching Covariate
region Factor (w/4 levels) Region of the hospital: Northeast, Midwest, South and West Covariate
ownership Factor (w/3 levels) Ownership of the hospital: Government, Private, non-profit and Private for-profit Covariate
year Character Year of the procedure Covariate
vol_reduction numeric Type of lung volume reduction: bronchoscopic (1) or surgical (0) Intervention/Comparator
vol_reduction_f Factor (w/2 levels) Type of lung volume reduction: bronchoscopic BLVR or surgical LVRS Intervention/Comparator
arf Numeric Presence of acute respiratory failure: 0/1 Outcome 1
arf_f Factor (w/2 levels) Presence of acute respiratory failure: Yes or No Outcome 1
ptx Numeric Presence of pneumothorax: 0/1 Outcome 1
ptx_f Factor (w/2 levels) Presence of pneumothorax: Yes or No Outcome 2

11.2 Creating Table (1)

Code
blvr.vars <- c("age", "race", "sex", "crf", "income", "insurance",
              "hospital_size", "teaching", "region", "ownership", "year")
              
blvr.reduction <- c("vol_reduction_f")

blvr_table1 <- CreateTableOne(data = blvr, 
                       vars = blvr.vars,
                       strata = blvr.reduction,
                       test = FALSE)
Code
print(blvr_table1) |> kable() |> kable_material()
                             Stratified by vol_reduction_f
                              LVRS          BLVR         
  n                             265           253        
  age (mean (SD))             56.40 (15.07) 67.88 (8.39) 
  race (%)                                               
     White                      204 (77.0)    228 (90.1) 
     Black                       31 (11.7)     17 ( 6.7) 
     Other                       30 (11.3)      8 ( 3.2) 
  sex = Female (%)               76 (28.7)    110 (43.5) 
  crf = Yes (%)                  34 (12.8)    117 (46.2) 
  income (%)                                             
     0-25th percentile           74 (27.9)     68 (26.9) 
     26th to 50th percentile     68 (25.7)     68 (26.9) 
     51st to 75th percentile     70 (26.4)     61 (24.1) 
     76th to 100th percentile    53 (20.0)     56 (22.1) 
  insurance (%)                                          
     Medicare                   108 (40.8)    196 (77.5) 
     Medicaid                    51 (19.2)     11 ( 4.3) 
     Private                    106 (40.0)     46 (18.2) 
  hospital_size (%)                                      
     Small                       12 ( 4.5)     27 (10.7) 
     Medium                      55 (20.8)     42 (16.6) 
     Large                      198 (74.7)    184 (72.7) 
  teaching (%)                                           
     Rural                        7 ( 2.6)      7 ( 2.8) 
     Urban nonteaching           20 ( 7.5)     28 (11.1) 
     Urban teaching             238 (89.8)    218 (86.2) 
  region (%)                                             
     Northeast                   63 (23.8)     32 (12.6) 
     Midwest                     79 (29.8)     77 (30.4) 
     South                       78 (29.4)    105 (41.5) 
     West                        45 (17.0)     39 (15.4) 
  ownership (%)                                          
     Government                  37 (14.0)     31 (12.3) 
     Private, non-profit        210 (79.2)    214 (84.6) 
     Private, for-profit         18 ( 6.8)      8 ( 3.2) 
  year (%)                                               
     2018                        85 (32.1)      6 ( 2.4) 
     2019                        97 (36.6)     76 (30.0) 
     2020                        83 (31.3)    171 (67.6) 
LVRS BLVR
n 265 253
age (mean (SD)) 56.40 (15.07) 67.88 (8.39)
race (%)
White 204 (77.0) 228 (90.1)
Black 31 (11.7) 17 ( 6.7)
Other 30 (11.3) 8 ( 3.2)
sex = Female (%) 76 (28.7) 110 (43.5)
crf = Yes (%) 34 (12.8) 117 (46.2)
income (%)
0-25th percentile 74 (27.9) 68 (26.9)
26th to 50th percentile 68 (25.7) 68 (26.9)
51st to 75th percentile 70 (26.4) 61 (24.1)
76th to 100th percentile 53 (20.0) 56 (22.1)
insurance (%)
Medicare 108 (40.8) 196 (77.5)
Medicaid 51 (19.2) 11 ( 4.3)
Private 106 (40.0) 46 (18.2)
hospital_size (%)
Small 12 ( 4.5) 27 (10.7)
Medium 55 (20.8) 42 (16.6)
Large 198 (74.7) 184 (72.7)
teaching (%)
Rural 7 ( 2.6) 7 ( 2.8)
Urban nonteaching 20 ( 7.5) 28 (11.1)
Urban teaching 238 (89.8) 218 (86.2)
region (%)
Northeast 63 (23.8) 32 (12.6)
Midwest 79 (29.8) 77 (30.4)
South 78 (29.4) 105 (41.5)
West 45 (17.0) 39 (15.4)
ownership (%)
Government 37 (14.0) 31 (12.3)
Private, non-profit 210 (79.2) 214 (84.6)
Private, for-profit 18 ( 6.8) 8 ( 3.2)
year (%)
2018 85 (32.1) 6 ( 2.4)
2019 97 (36.6) 76 (30.0)
2020 83 (31.3) 171 (67.6)

Based on table 1, patients receiving bronchoscopic lung volume reduction (BLVR) were noted to be older (mean age of 67 years), with a higher percentage of females (43.5%), likely due to a preference of a less invasive approach compared to the surgical group. BLVR patients also tended to have a higher percentage of chronic respiratory failure, likely due to the ease of valve removal should the patient’s baseline respiratory failure worsen after procedure. The patient population was more white (90%), with medicare as the primary insurance provider (77.5%) which is likely due to the age bracket in this group (65 years and older).

As expected both groups had a similar distribution of median income per zip code, and majority of these procedures were done in large, non for profit and urban teaching hospitals. Lung volume reduction surgery number distribution was steady over the years included (roughly around 30%), however bronchoscopic lung volume reduction numbers significantly increased over the years (starting in 2019 and 2020)

12 Unadjusted Outcomes

12.1 Acute Respiratory Failure arf

Code
ggplot(blvr, aes(x = vol_reduction_f, fill = arf_f)) +
    geom_bar(aes(y = (after_stat(count))/sum(after_stat(count))), position = "dodge") +
    scale_y_continuous(labels = scales::percent) +
    scale_fill_viridis_d(option = "turbo", alpha = 0.5) +  labs(title = "Respiratory Failure", x = "Volume Reduction", y = "% of Patients", fill = "Resp. Failure")

Code
blvr |> tabyl(arf_f, vol_reduction_f) |> kable() |> kable_material()
arf_f LVRS BLVR
No 235 209
Yes 30 44

The visualization above suggested that BLVR patients tended to suffer a higher percentage of acute respiratory failure

A logistic regression model was fit to assess the odds of respiratory failure based on the type of volume reduction received, then the coefficient was exponentiated to obtain an odds ratio

Code
unadjusted_arf <- glm(arf ~ vol_reduction, data=blvr, family=binomial)
Code
kable(tidy(unadjusted_arf, conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE), digits = 2) |> kable_paper()
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 0.13 0.19 -10.62 0.00 0.09 0.18
vol_reduction 1.65 0.26 1.96 0.05 1.00 2.74

If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR has a 1.65 odds of acute respiratory failure (95% CI (1,2.74)). This confidence interval contained 1, thus it was difficult to reliably say that BLVR would meaningfully lead to a higher odds of acute respiratory failure

Overall, it seems that the patients receiving BLVR tended to have higher odds of respiratory failure, even though the odds ratio’s interval contained 1, this odds ratio shows an interesting trend

12.2 Pneumothorax ptx

Code
ggplot(blvr, aes(x = vol_reduction_f, fill = ptx_f)) +
    geom_bar(aes(y = (after_stat(count))/sum(after_stat(count))), position = "dodge") +
    scale_y_continuous(labels = scales::percent) +
    scale_fill_viridis_d(option = "turbo", alpha = 0.5) +
  labs(title = "Pneumothorax", x = "Volume Reduction", y = "% of Patients", fill = "PTX")

Code
blvr |> tabyl(ptx_f, vol_reduction_f) |> kable() |> kable_material()
ptx_f LVRS BLVR
No 159 200
Yes 106 53

The visualization above suggested that BLVR patients tended to suffer a lower percentage of pneumothorax

A logistic regression model was fit to assess the odds of pneumothorax based on the type of volume reduction received, then the coefficient was exponentiated to obtain an odds ratio

Code
unadjusted_ptx <- glm(ptx ~ vol_reduction, data=blvr, family=binomial)
Code
kable(tidy(unadjusted_ptx, conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE), digits = 2) |> kable_paper()
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 0.67 0.13 -3.23 0 0.52 0.85
vol_reduction 0.40 0.20 -4.64 0 0.27 0.58

If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR had a 0.4 odds of pneumothorax (95% CI (0.27,0.58)), which is lower odds than the subject receiving LVRS, which is plausible as LVRS is more invasive.

13 Propensity Score

13.1 Fitting the model

A logistic regression model was used to estimate the propensity score. Having BLVR/LVRS was modeled as the outcome, with patient and hospital characteristics as independent variables. Patient characteristics used in the propensity score included age, sex, race, crf, income and insurance. Hospital characteristics included hospital_size, teaching status, region, ownership and year of the procedure. The raw and linear propensity score values were saved to the blvr tibble.

Code
psmodel <- glm(vol_reduction ~ age + sex + race + crf + income + insurance + 
                 hospital_size + teaching + region + ownership + year,
               family=binomial, data=blvr)
summary(psmodel)

Call:
glm(formula = vol_reduction ~ age + sex + race + crf + income + 
    insurance + hospital_size + teaching + region + ownership + 
    year, family = binomial, data = blvr)

Coefficients:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    -4.696344   1.503545  -3.124 0.001787 ** 
age                             0.064113   0.015710   4.081 4.48e-05 ***
sexFemale                       0.716995   0.258418   2.775 0.005528 ** 
raceBlack                      -0.481109   0.451914  -1.065 0.287055    
raceOther                      -1.521225   0.599723  -2.537 0.011195 *  
crfYes                          1.294080   0.288002   4.493 7.01e-06 ***
income26th to 50th percentile   0.393352   0.350994   1.121 0.262423    
income51st to 75th percentile  -0.286941   0.342652  -0.837 0.402361    
income76th to 100th percentile  0.026792   0.386403   0.069 0.944722    
insuranceMedicaid              -1.353403   0.515379  -2.626 0.008639 ** 
insurancePrivate               -0.469448   0.318487  -1.474 0.140484    
hospital_sizeMedium            -2.400828   0.671416  -3.576 0.000349 ***
hospital_sizeLarge             -2.077553   0.616456  -3.370 0.000751 ***
teachingUrban nonteaching       0.678464   0.954277   0.711 0.477102    
teachingUrban teaching         -0.472681   0.823311  -0.574 0.565885    
regionMidwest                   0.044408   0.388618   0.114 0.909021    
regionSouth                     0.668451   0.378964   1.764 0.077750 .  
regionWest                      0.167103   0.437435   0.382 0.702457    
ownershipPrivate, non-profit   -0.006272   0.357697  -0.018 0.986011    
ownershipPrivate, for-profit   -1.513048   0.774833  -1.953 0.050850 .  
year2019                        2.211839   0.517582   4.273 1.93e-05 ***
year2020                        3.281191   0.507832   6.461 1.04e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 717.82  on 517  degrees of freedom
Residual deviance: 436.37  on 496  degrees of freedom
AIC: 480.37

Number of Fisher Scoring iterations: 6
Code
blvr$ps <- psmodel$fitted
blvr$linps <- psmodel$linear.predictors

13.2 Describing the overlap

Code
ggplot(blvr, aes(x = vol_reduction_f, y = ps, color = vol_reduction_f)) + 
    geom_boxplot() +
    geom_jitter(width = 0.1, alpha = 0.3) +
    scale_color_viridis_d(option = "plasma", end = 0.5) + 
  labs(title = "Boxplot & Jitter", subtitle = "Distribution Of the Raw PS",  x = "Type of Procedure", y = "Propensity Scores",
       caption = "* Line indicates sample median", col = "Procedure")

Code
ggplot(blvr, aes(x = linps, fill = vol_reduction_f)) +
    geom_density(alpha = 0.3) + 
  scale_fill_viridis_d(option = "turbo") +
  labs(title = "Density Plot", subtitle = "For the Linear Propensity Score",  x = "Linps", y = "Density",
       fill = "Volume Reduction")

Code
blvr |> group_by(vol_reduction_f) |> skim_without_charts(ps, linps)
Data summary
Name group_by(blvr, vol_reduct…
Number of rows 518
Number of columns 20
_______________________
Column type frequency:
numeric 2
________________________
Group variables vol_reduction_f

Variable type: numeric

skim_variable vol_reduction_f n_missing complete_rate mean sd p0 p25 p50 p75 p100
ps LVRS 0 1 0.27 0.26 0.00 0.03 0.18 0.47 0.93
ps BLVR 0 1 0.72 0.24 0.03 0.60 0.80 0.91 1.00
linps LVRS 0 1 -1.81 2.14 -7.45 -3.37 -1.53 -0.13 2.52
linps BLVR 0 1 1.29 1.53 -3.61 0.40 1.38 2.27 5.64

The overlap between the propensity scores in intervention groups was assessed numerically and graphically. The boxplot and the density plots show some overlap in the propensity scores between the intervention and the treatment.

The LVRS group had a lower propensity score mean, with problematic values falling below 0.01. The BLVR group also had some problematic high values over 0.99. Troubleshooting will be addressed at a later step

13.3 Unadjusted Rubin’s Rules

Code
unadjusted_rubin1 <- 
  abs(100*(mean(blvr$linps[blvr$vol_reduction==1]) -
             mean(blvr$linps[blvr$vol_reduction==0])) / 
                           sd(blvr$linps))

unadjusted_rubin2 <- with(blvr, var(linps[vol_reduction == 1]) / var(linps[vol_reduction == 0]))
Code
unadjusted_rubin1; unadjusted_rubin2
[1] 127.7055
[1] 0.506223

Rubin’s rule 1 value was higher than 50%, this meant that running an unadjusted model could not be justified without applying the propensity scores to account for selection bias.

Rubin’s rule 2 was on the lower bound of acceptable value as well, as the variance ratio of the linear propensity score in the BLVR (intervention) group to the variance of the linear propensity score in the LVRS (control) group should be close to 1 (and ideally between 0.8 and 1.25)

13.4 Troubleshooting the Extreme Propensity Score Values

Code
blvr |> count(ps >= 0.99); blvr |> count(ps < 0.01)
# A tibble: 2 × 2
  `ps >= 0.99`     n
  <lgl>        <int>
1 FALSE          514
2 TRUE             4
# A tibble: 2 × 2
  `ps < 0.01`     n
  <lgl>       <int>
1 FALSE         490
2 TRUE           28

There were twenty eight observations with propensity score values below 0.01 and four observations with propensity score values above 0.99. The extreme lows seemed to be present in the LVRS group while the extreme highs were present in the BLVR group

Code
mosaic::favstats(ps ~ vol_reduction_f, data = blvr)
Registered S3 method overwritten by 'mosaic':
  method                           from   
  fortify.SpatialPolygonsDataFrame ggplot2
  vol_reduction_f          min         Q1    median        Q3       max
1            LVRS 0.0005808267 0.03308973 0.1778799 0.4680281 0.9256507
2            BLVR 0.0263071482 0.59915004 0.7983214 0.9061054 0.9964521
       mean        sd   n missing
1 0.2679640 0.2584781 265       0
2 0.7193262 0.2402497 253       0

The following strategies were attempted in order to resolve the above issue:

  • Colinearity check using vif didn’t reveal problematic colinearity between the variables used in the propensity score model
Code
car::vif(psmodel)
                  GVIF Df GVIF^(1/(2*Df))
age           1.634643  1        1.278532
sex           1.070972  1        1.034878
race          1.204488  2        1.047612
crf           1.105123  1        1.051248
income        1.428943  3        1.061294
insurance     1.667039  2        1.136283
hospital_size 1.403948  2        1.088523
teaching      1.250028  2        1.057377
region        1.470981  3        1.066435
ownership     1.355067  2        1.078922
year          1.164515  2        1.038810
  • Re-constructing the propensity score model, adding one covariate at a time, while checking which variable drives the propensity score values below 0.1 or above 0.99 showed that each of the following covariates race, year, sex, ownership and crf could independently drop the propensity score below 0.01 or drive it up above 0.99. Excluding these covariates form the model would be problematic in itself

  • Collapsing multicategorical variables was attempted, these variables included year, race, ownership, insurance, hospital_size, teaching and region. It did not result in the improvement of the propensity score even when each covariate was added individually after collapsing.

Filtering obaservations with problematic propensity scores proved to be the most reasonable approach in this situation. This resulted in a sample with 237 observations on the control intervention (LVRS) and 249 observation on the studied intervention (BLVR). The resulting dataset blvr_matching was used for propensity score matching analyses

Code
blvr_matching <- blvr |> filter(ps <= 0.990 & ps >= 0.010)
Code
mosaic::favstats(ps ~ vol_reduction_f, data = blvr_matching)
  vol_reduction_f        min         Q1    median        Q3       max      mean
1            LVRS 0.01020130 0.05651206 0.2244966 0.4944178 0.9256507 0.2991748
2            BLVR 0.02630715 0.58669227 0.7927773 0.9014236 0.9892394 0.7149098
         sd   n missing
1 0.2558878 237       0
2 0.2396084 249       0

A density plot was created to describe the overlap of the linear propensity scores between the two intervention groups after balancing the dataset, the overlap is described graphically in the density plot below.

Code
ggplot(blvr_matching, aes(x = linps, fill = vol_reduction_f)) +
    geom_density(alpha = 0.3) + 
  scale_fill_viridis_d(option = "turbo") +
  labs(title = "Density Plot", subtitle = "For the Linear Propensity Score",  x = "Linps", y = "Density",
       fill = "Volume Reduction")

14 Approach (1): Matching

Propensity score matching was conducted using a 1:1 caliper matching on the linear propensity score without replacement from the Matching package. A default caliper value of 0.2 was selected. The type of matching was chosen based on the similar sample size between the intervention and control groups.

Code
X <- blvr_matching$linps
Tr <- as.logical(blvr_matching$vol_reduction)
match <- Match(Tr = Tr, X = X, M = 1, 
                caliper = 0.2,
                estimand = "ATT", replace = FALSE, 
                ties = FALSE)
Warning in Match(Tr = Tr, X = X, M = 1, caliper = 0.2, estimand = "ATT", :
replace==FALSE, but there are more (weighted) treated obs than control obs.
Some treated obs will not be matched.  You may want to estimate ATC instead.
Code
summary(match)

Estimate...  0 
SE.........  0 
T-stat.....  NaN 
p.val......  NA 

Original number of observations..............  486 
Original number of treated obs...............  249 
Matched number of observations...............  107 
Matched number of observations  (unweighted).  107 

Caliper (SDs)........................................   0.2 
Number of obs dropped by 'exact' or 'caliper'  142 

One hundred and seven BLVR patients were matched to one hundred and seven LVRS patients. The degree of covariate imbalance was assessed before and after propensity score matching for each covariate as well as the linear and raw propensity scores

Code
covs <- blvr_matching |> 
  select(age, sex, race, crf, insurance, income, hospital_size, region, teaching, ownership, year, ps, linps)

b <- bal.tab(match,
             treat = blvr_matching$vol_reduction,
             covs = covs,
             stats = c("m", "v"),
             quick = FALSE, un = TRUE, disp.v.ratio = TRUE)
b
Balance Measures
                                   Type Diff.Un V.Ratio.Un Diff.Adj V.Ratio.Adj
age                             Contin.  1.1098     0.3830   0.1295      0.8005
sex_Female                       Binary  0.1135          .  -0.0280           .
race_White                       Binary  0.0977          .   0.0093           .
race_Black                       Binary -0.0497          .   0.0093           .
race_Other                       Binary -0.0480          .  -0.0187           .
crf_Yes                          Binary  0.3104          .   0.0187           .
insurance_Medicare               Binary  0.3238          .   0.0841           .
insurance_Medicaid               Binary -0.1204          .  -0.0093           .
insurance_Private                Binary -0.2034          .  -0.0748           .
income_0-25th percentile         Binary -0.0052          .   0.0000           .
income_26th to 50th percentile   Binary -0.0008          .  -0.0280           .
income_51st to 75th percentile   Binary -0.0333          .   0.0280           .
income_76th to 100th percentile  Binary  0.0392          .   0.0000           .
hospital_size_Small              Binary  0.0542          .  -0.0187           .
hospital_size_Medium             Binary -0.0170          .   0.0000           .
hospital_size_Large              Binary -0.0372          .   0.0187           .
region_Northeast                 Binary -0.0951          .   0.0374           .
region_Midwest                   Binary  0.0181          .   0.0187           .
region_South                     Binary  0.0934          .   0.0000           .
region_West                      Binary -0.0164          .  -0.0561           .
teaching_Rural                   Binary -0.0014          .   0.0093           .
teaching_Urban nonteaching       Binary  0.0367          .  -0.0093           .
teaching_Urban teaching          Binary -0.0353          .   0.0000           .
ownership_Government             Binary -0.0314          .  -0.0467           .
ownership_Private, non-profit    Binary  0.0499          .   0.0280           .
ownership_Private, for-profit    Binary -0.0185          .   0.0187           .
year_2018                        Binary -0.2333          .  -0.0467           .
year_2019                        Binary -0.0956          .   0.1402           .
year_2020                        Binary  0.3289          .  -0.0935           .
ps                              Contin.  1.7351     0.8768   0.0429      1.0478
linps                           Contin.  1.7602     0.7249   0.0326      1.0488

Sample sizes
          Control Treated
All           237     249
Matched       107     107
Unmatched     130       0
Discarded       0     142
Code
bal.plot(match,
             treat = blvr_matching$vol_reduction,
             covs = covs,
             var.name = "linps",
             which = "both",
             sample.names = 
               c("Unmatched Sample", "Matched Sample")) +
  scale_fill_viridis_d(option = "turbo") + 
  labs(title = "Distributional Balance for linear PS", subtitle = "1:1 Caliper Matching", x = "Linear PS") 

A Love plot of standardized differences before and after matching was obtained to evaluate matching success, showing improvement over the unmatched sample

Code
love.plot(b, threshold = .1, size = 3,
             var.order = "unadjusted",
             colors = c("hotpink3", "purple4"),
             stars = "raw") + 
  labs(title = "Love Plot of Standardized Differences", subtitle = "1:1 Caliper Matching")

Rubin’s first and second rules were assessed after matching, showing improvement over the unadjusted Rubin’s rules.

Code
matches <- factor(rep(match$index.treated, 2))

blvr.matchedsample <- cbind(matches, blvr_matching[c(match$index.control, match$index.treated),])
Code
matched_rubin1 <- with(blvr.matchedsample,
                      abs(100*(mean(linps[vol_reduction==1])-mean(linps[vol_reduction==0]))/sd(linps)))

matched_rubin2 <- with(blvr.matchedsample, var(linps[vol_reduction==1])/var(linps[vol_reduction==0]))
Code
matched_rubin1; matched_rubin2
[1] 4.060997
[1] 1.048778

Rubin’s rule 1 was less than 10%, and Rubin’s rule 2 was around 1 which was ideal.

15 Outcomes after Matching

15.1 Acute respiratory failure arf

Code
matched_arf <- clogit(arf ~ vol_reduction + strata(matches), data=blvr.matchedsample)
Code
kable(tidy(matched_arf, conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE), digits = 2) |> kable_material()
term estimate std.error statistic p.value conf.low conf.high
vol_reduction 1.73 0.38 1.44 0.15 0.82 3.63

If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR has a 1.73 odds of acute respiratory failure (95% CI (0.82,3.63)). The subject receiving BLVR tended to have a higher odds ratio of respiratory failue in comparison to patients receiving LVRS, however this confidence interval crosses 1

15.2 Pneumothorax ptx

Code
matched_ptx <- clogit(ptx ~ vol_reduction + strata(matches), data=blvr.matchedsample)
Code
kable(tidy(matched_ptx, exponentiate = TRUE, conf.level = 0.95, conf.int = TRUE), digits = 2) |> kable_material()
term estimate std.error statistic p.value conf.low conf.high
vol_reduction 0.68 0.3 -1.3 0.19 0.38 1.22

If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR has a 0.68 odds of pneumothorax (95% CI (0.38,1.22)). Patients receiving BLVR seemed to have a higher odds ratio of pneumothorax in comparison to patients receiving BLVR, however this confidence interval contained 1, the results were not statistically detectable and a stability analysis instead of a sensitivity analysis was carried out.

15.2.1 Stability Analysis

  • Trying different imputed datasets (generated through mice package) resulted in similar outcome estimates (odds ratio)

  • The matching strategies described below were attempted. None of these strategies provided an improvement in the effective sample size, propensity score distributional balance, standardized differences or Rubin’s rules in comparison to the matching strategy used in this analysis (1:1 matching with a caliper)

Match Trtmnt Ctrl Distributional Balance Std Differences Rubin’s 1 Rubin’s 2
1:1 Greedy W replacement 249 75 Better overlap, not ideal better but dispersed 8.8 1.27
Genetic matching 273 78 Better overlap, not ideal Dispersed 8.8 1.26
1:1 Neartest neighbor matching W Replacement 249 79 Better overlap, not ideal better but dispersed (> 0.1) 8 1.21
  • The following steps attempted comparing the odds ratios obtained through matching to the odds ratios after ATT weighting and a double robust approach (by adjusting for the linear propensity score in the regression model)

16 Approach (2): Propensity Score Weighting

Given the initial difference in sample sizes across the two volume reduction cohorts, propensity score matching has left a selection of the patients unmatched. Weighting by the inverse propensity score (average of the treatment effect on the treated) was next carried out using the initial database blvr (the database that still contained propensity scores below 0.01 and above 0.99)

Code
blvr$wts1 <- ifelse(blvr$vol_reduction==1, 1, blvr$ps/(1-blvr$ps))

Then the resulting ATT (average treatment effect on the treated) weights were plotted in the figure below, where the treatment group got the same weight of 1, and the control group had varying weights

Code
ggplot(blvr, aes(x = ps, y = wts1, color = vol_reduction_f)) +
    geom_point(size = 2, alpha = 0.4) + 
    facet_wrap(~ vol_reduction_f) +
    labs(x = "Estimated Propensity for Volume Reduction",
         y = "ATT weights for blvr",
         title = "ATT weighting structure",
         subtitle = "In the BLVR dataset",
         col = "Volume Reduction")

The balance after weighting was assessed using dx.wts from the twang package. Prior to weighting, the control sample was 237 participants, and after weighting, the effective sample size was 54 participants, representing a sample size that would be required to achieve the same level of precision if that sample were a simple random sample

Code
blvr_df <- data.frame(blvr)

covlist <- c("age", "sex", "race", "crf", "insurance", "income", "hospital_size", "teaching", "region", "ownership", "year", "ps", "linps") 
  
bal.wts1 <- dx.wts(x=blvr_df$wts1, data=blvr_df, vars=covlist, 
                   treat.var="vol_reduction", estimand="ATT")
bal.wts1
  type n.treat n.ctrl ess.treat ess.ctrl    max.es   mean.es    max.ks
1  unw     253    265       253 265.0000 2.0306797 0.4589984 0.6194496
2          253    265       253  54.7126 0.3077257 0.1118164 0.2134387
     mean.ks iter
1 0.14873095   NA
2 0.05053566   NA

The standardized differences were compared before and after weighting using a Love.plot. The weighting process showed a more satisfactory balance of standardized differences in comparison to the unweighted standardized differences.

Code
bal.table(bal.wts1)
$unw
                                 tx.mn tx.sd  ct.mn  ct.sd std.eff.sz   stat
age                             67.881 8.390 56.400 15.065      1.368 10.789
sex:Male                         0.565 0.496  0.713  0.452     -0.299 12.293
sex:Female                       0.435 0.496  0.287  0.452      0.299     NA
race:White                       0.901 0.298  0.770  0.421      0.440  8.925
race:Black                       0.067 0.250  0.117  0.321     -0.199     NA
race:Other                       0.032 0.175  0.113  0.317     -0.466     NA
crf:No                           0.538 0.499  0.872  0.334     -0.670 69.838
crf:Yes                          0.462 0.499  0.128  0.334      0.670     NA
insurance:Medicare               0.775 0.418  0.408  0.491      0.879 37.291
insurance:Medicaid               0.043 0.204  0.192  0.394     -0.731     NA
insurance:Private                0.182 0.386  0.400  0.490     -0.566     NA
income:0-25th percentile         0.269 0.443  0.279  0.449     -0.024  0.225
income:26th to 50th percentile   0.269 0.443  0.257  0.437      0.027     NA
income:51st to 75th percentile   0.241 0.428  0.264  0.441     -0.054     NA
income:76th to 100th percentile  0.221 0.415  0.200  0.400      0.051     NA
hospital_size:Small              0.107 0.309  0.045  0.208      0.199  3.868
hospital_size:Medium             0.166 0.372  0.208  0.406     -0.112     NA
hospital_size:Large              0.727 0.445  0.747  0.435     -0.045     NA
teaching:Rural                   0.028 0.164  0.026  0.160      0.008  0.965
teaching:Urban nonteaching       0.111 0.314  0.075  0.264      0.112     NA
teaching:Urban teaching          0.862 0.345  0.898  0.302     -0.106     NA
region:Northeast                 0.126 0.332  0.238  0.426     -0.335  4.752
region:Midwest                   0.304 0.460  0.298  0.457      0.014     NA
region:South                     0.415 0.493  0.294  0.456      0.245     NA
region:West                      0.154 0.361  0.170  0.375     -0.043     NA
ownership:Government             0.123 0.328  0.140  0.347     -0.052  2.065
ownership:Private, non-profit    0.846 0.361  0.792  0.406      0.148     NA
ownership:Private, for-profit    0.032 0.175  0.068  0.252     -0.207     NA
year:2018                        0.024 0.152  0.321  0.467     -1.952 50.600
year:2019                        0.300 0.458  0.366  0.482     -0.143     NA
year:2020                        0.676 0.468  0.313  0.464      0.775     NA
ps                               0.719 0.240  0.268  0.258      1.879 20.616
linps                            1.288 1.526 -1.811  2.145      2.031 19.031
                                    p    ks ks.pval
age                             0.000 0.359   0.000
sex:Male                        0.000 0.148   0.000
sex:Female                         NA 0.148   0.000
race:White                      0.000 0.131   0.000
race:Black                         NA 0.050   0.000
race:Other                         NA 0.082   0.000
crf:No                          0.000 0.334   0.000
crf:Yes                            NA 0.334   0.000
insurance:Medicare              0.000 0.367   0.000
insurance:Medicaid                 NA 0.149   0.000
insurance:Private                  NA 0.218   0.000
income:0-25th percentile        0.879 0.010   0.879
income:26th to 50th percentile     NA 0.012   0.879
income:51st to 75th percentile     NA 0.023   0.879
income:76th to 100th percentile    NA 0.021   0.879
hospital_size:Small             0.021 0.061   0.021
hospital_size:Medium               NA 0.042   0.021
hospital_size:Large                NA 0.020   0.021
teaching:Rural                  0.381 0.001   0.381
teaching:Urban nonteaching         NA 0.035   0.381
teaching:Urban teaching            NA 0.036   0.381
region:Northeast                0.003 0.111   0.003
region:Midwest                     NA 0.006   0.003
region:South                       NA 0.121   0.003
region:West                        NA 0.016   0.003
ownership:Government            0.127 0.017   0.127
ownership:Private, non-profit      NA 0.053   0.127
ownership:Private, for-profit      NA 0.036   0.127
year:2018                       0.000 0.297   0.000
year:2019                          NA 0.066   0.000
year:2020                          NA 0.363   0.000
ps                              0.000 0.619   0.000
linps                           0.000 0.619   0.000

[[2]]
                                 tx.mn tx.sd  ct.mn ct.sd std.eff.sz  stat
age                             67.881 8.390 66.767 8.471      0.133 0.923
sex:Male                         0.565 0.496  0.613 0.487     -0.096 0.372
sex:Female                       0.435 0.496  0.387 0.487      0.096    NA
race:White                       0.901 0.298  0.866 0.341      0.118 0.576
race:Black                       0.067 0.250  0.104 0.306     -0.148    NA
race:Other                       0.032 0.175  0.030 0.170      0.010    NA
crf:No                           0.538 0.499  0.655 0.475     -0.236 2.117
crf:Yes                          0.462 0.499  0.345 0.475      0.236    NA
insurance:Medicare               0.775 0.418  0.782 0.413     -0.018 0.111
insurance:Medicaid               0.043 0.204  0.034 0.182      0.046    NA
insurance:Private                0.182 0.386  0.184 0.387     -0.005    NA
income:0-25th percentile         0.269 0.443  0.325 0.469     -0.128 0.790
income:26th to 50th percentile   0.269 0.443  0.180 0.384      0.200    NA
income:51st to 75th percentile   0.241 0.428  0.265 0.442     -0.057    NA
income:76th to 100th percentile  0.221 0.415  0.229 0.420     -0.019    NA
hospital_size:Small              0.107 0.309  0.101 0.301      0.019 0.006
hospital_size:Medium             0.166 0.372  0.169 0.375     -0.009    NA
hospital_size:Large              0.727 0.445  0.730 0.444     -0.006    NA
teaching:Rural                   0.028 0.164  0.017 0.130      0.063 3.263
teaching:Urban nonteaching       0.111 0.314  0.037 0.188      0.236    NA
teaching:Urban teaching          0.862 0.345  0.946 0.226     -0.245    NA
region:Northeast                 0.126 0.332  0.208 0.406     -0.247 1.450
region:Midwest                   0.304 0.460  0.334 0.472     -0.065    NA
region:South                     0.415 0.493  0.363 0.481      0.105    NA
region:West                      0.154 0.361  0.094 0.292      0.166    NA
ownership:Government             0.123 0.328  0.099 0.299      0.071 1.134
ownership:Private, non-profit    0.846 0.361  0.887 0.316     -0.115    NA
ownership:Private, for-profit    0.032 0.175  0.013 0.115      0.104    NA
year:2018                        0.024 0.152  0.032 0.176     -0.055 0.187
year:2019                        0.300 0.458  0.318 0.466     -0.039    NA
year:2020                        0.676 0.468  0.650 0.477      0.056    NA
ps                               0.719 0.240  0.663 0.223      0.233 1.767
linps                            1.288 1.526  0.818 1.184      0.308 2.489
                                    p    ks ks.pval
age                             0.356 0.099   0.734
sex:Male                        0.542 0.048   0.542
sex:Female                         NA 0.048   0.542
race:White                      0.552 0.035   0.552
race:Black                         NA 0.037   0.552
race:Other                         NA 0.002   0.552
crf:No                          0.146 0.118   0.146
crf:Yes                            NA 0.118   0.146
insurance:Medicare              0.870 0.008   0.870
insurance:Medicaid                 NA 0.009   0.870
insurance:Private                  NA 0.002   0.870
income:0-25th percentile        0.494 0.057   0.494
income:26th to 50th percentile     NA 0.089   0.494
income:51st to 75th percentile     NA 0.024   0.494
income:76th to 100th percentile    NA 0.008   0.494
hospital_size:Small             0.990 0.006   0.990
hospital_size:Medium               NA 0.003   0.990
hospital_size:Large                NA 0.003   0.990
teaching:Rural                  0.044 0.010   0.044
teaching:Urban nonteaching         NA 0.074   0.044
teaching:Urban teaching            NA 0.084   0.044
region:Northeast                0.229 0.082   0.229
region:Midwest                     NA 0.030   0.229
region:South                       NA 0.052   0.229
region:West                        NA 0.060   0.229
ownership:Government            0.319 0.023   0.319
ownership:Private, non-profit      NA 0.042   0.319
ownership:Private, for-profit      NA 0.018   0.319
year:2018                       0.784 0.008   0.784
year:2019                          NA 0.018   0.784
year:2020                          NA 0.026   0.784
ps                              0.078 0.213   0.029
linps                           0.013 0.213   0.029
Code
bal.before.wts1 <- bal.table(bal.wts1)[1]
bal.after.wts1 <- bal.table(bal.wts1)[2]

balance.att.weights <- data.frame(names = rownames(bal.before.wts1$unw), 
                              pre.weighting = 100*bal.before.wts1$unw$std.eff.sz, 
                              ATT.weighted = 100*bal.after.wts1[[1]]$std.eff.sz)

balance.att.weights <- gather(balance.att.weights, timing, szd, 2:3)

balance.att.weights
                             names        timing    szd
1                              age pre.weighting  136.8
2                         sex:Male pre.weighting  -29.9
3                       sex:Female pre.weighting   29.9
4                       race:White pre.weighting   44.0
5                       race:Black pre.weighting  -19.9
6                       race:Other pre.weighting  -46.6
7                           crf:No pre.weighting  -67.0
8                          crf:Yes pre.weighting   67.0
9               insurance:Medicare pre.weighting   87.9
10              insurance:Medicaid pre.weighting  -73.1
11               insurance:Private pre.weighting  -56.6
12        income:0-25th percentile pre.weighting   -2.4
13  income:26th to 50th percentile pre.weighting    2.7
14  income:51st to 75th percentile pre.weighting   -5.4
15 income:76th to 100th percentile pre.weighting    5.1
16             hospital_size:Small pre.weighting   19.9
17            hospital_size:Medium pre.weighting  -11.2
18             hospital_size:Large pre.weighting   -4.5
19                  teaching:Rural pre.weighting    0.8
20      teaching:Urban nonteaching pre.weighting   11.2
21         teaching:Urban teaching pre.weighting  -10.6
22                region:Northeast pre.weighting  -33.5
23                  region:Midwest pre.weighting    1.4
24                    region:South pre.weighting   24.5
25                     region:West pre.weighting   -4.3
26            ownership:Government pre.weighting   -5.2
27   ownership:Private, non-profit pre.weighting   14.8
28   ownership:Private, for-profit pre.weighting  -20.7
29                       year:2018 pre.weighting -195.2
30                       year:2019 pre.weighting  -14.3
31                       year:2020 pre.weighting   77.5
32                              ps pre.weighting  187.9
33                           linps pre.weighting  203.1
34                             age  ATT.weighted   13.3
35                        sex:Male  ATT.weighted   -9.6
36                      sex:Female  ATT.weighted    9.6
37                      race:White  ATT.weighted   11.8
38                      race:Black  ATT.weighted  -14.8
39                      race:Other  ATT.weighted    1.0
40                          crf:No  ATT.weighted  -23.6
41                         crf:Yes  ATT.weighted   23.6
42              insurance:Medicare  ATT.weighted   -1.8
43              insurance:Medicaid  ATT.weighted    4.6
44               insurance:Private  ATT.weighted   -0.5
45        income:0-25th percentile  ATT.weighted  -12.8
46  income:26th to 50th percentile  ATT.weighted   20.0
47  income:51st to 75th percentile  ATT.weighted   -5.7
48 income:76th to 100th percentile  ATT.weighted   -1.9
49             hospital_size:Small  ATT.weighted    1.9
50            hospital_size:Medium  ATT.weighted   -0.9
51             hospital_size:Large  ATT.weighted   -0.6
52                  teaching:Rural  ATT.weighted    6.3
53      teaching:Urban nonteaching  ATT.weighted   23.6
54         teaching:Urban teaching  ATT.weighted  -24.5
55                region:Northeast  ATT.weighted  -24.7
56                  region:Midwest  ATT.weighted   -6.5
57                    region:South  ATT.weighted   10.5
58                     region:West  ATT.weighted   16.6
59            ownership:Government  ATT.weighted    7.1
60   ownership:Private, non-profit  ATT.weighted  -11.5
61   ownership:Private, for-profit  ATT.weighted   10.4
62                       year:2018  ATT.weighted   -5.5
63                       year:2019  ATT.weighted   -3.9
64                       year:2020  ATT.weighted    5.6
65                              ps  ATT.weighted   23.3
66                           linps  ATT.weighted   30.8
Code
ggplot(balance.att.weights, aes(x = szd, y = reorder(names, szd), color = timing)) +
    geom_point(size = 2, alpha = 0.5) + 
    geom_vline(xintercept = 0) +
    geom_vline(xintercept = c(-10,10), linetype = "dashed", col = "blue") +
    labs(x = "Standardized Differences", y = "", 
         title = "Love Plot of Standarized Diffs in The BLVR Data",
         subtitle = "Before and After ATT Weighting") 

Rubin’s rules were obtained after weighting. Rubin’s rule 1 was obtained by using the standardized effect size after weighting for the linear propensity score (0.28). When multiplied by 100, the result was 28% indicating that Rubin’s rule 1 was satisfied (although not perfectly). Rubin’s rule 2 was calculated by dividing the square root of treatment group standard deviation by the square root of control group standard deviation, resulting in a value of 1.54. Rubin’s rule 2 was also satisfied, although less perfectly.

17 Outcomes after Weighting

17.1 Acute Respiratory Failure arf

Code
blvr.design <- svydesign(ids=~1, weights=~wts1, data=blvr)
Code
weighted_arf <- svyglm(arf ~ vol_reduction, design=blvr.design, family=quasibinomial())
Code
kable(tidy(weighted_arf, conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE), digits = 2) |> kable_material() 
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 0.13 0.32 -6.49 0.00 0.07 0.24
vol_reduction 1.68 0.36 1.43 0.15 0.83 3.40

The estimate of treatment effect on acute respiratory failure was obtained by a regression model using survey weighting and a quasi binomial design.

The weighted odds ratio estimate for the BLVR group was 1.68, with a 95% confidence interval (0.83,3.4). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having acute respiratory failure for the person receiving BLVR was higher and 1.68 the odds of respiratory failure compared to the person receiving LVRS. It should be noted that this confidence interveal crosses 1

17.2 Pneumothorax ptx

Code
weighted_ptx <- svyglm(ptx ~ vol_reduction, design=blvr.design, family=quasibinomial())
Code
kable(tidy(weighted_ptx, conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE), digits = 2) |> kable_material() 
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 0.73 0.28 -1.13 0.26 0.43 1.26
vol_reduction 0.36 0.32 -3.21 0.00 0.19 0.67

The weighted odds ratio estimate for the `BLVR group was 0.36, with a 95% confidence interval (0.19,0.67). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having pneumothorax for the person receiving BLVR was lower and was 0.36 the odds of pneumothorax of the person receiving LVRS. The confidence interval did not cross 1

18 Apprach (3): A Double-Robust Estimate

A logistic regression for the outcomes of acute respiratory failure and pneumothorax was fitted with a regression adjustment for the linear propensity score to obtain a doubly-robust estimate and account for residual imbalance after weighting

Code
robust.design <- svydesign(ids=~1, weights=~wts1, data=blvr)
Code
robust_arf <- svyglm(arf ~ vol_reduction + linps, design=robust.design, family=quasibinomial())
Code
kable(tidy(robust_arf, exponentiate = TRUE, conf.int = TRUE, conf.level = 0.95), digits = 2) |> kable_material()
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 0.15 0.30 -6.23 0.00 0.08 0.28
vol_reduction 1.89 0.36 1.76 0.08 0.93 3.87
linps 0.75 0.08 -3.58 0.00 0.64 0.88

The doubly-robust odds ratio estimate for the BLVR group was 1.89, with a 95% confidence interval (0.93,3.87). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having acute respiratory failure for the person receiving BLVR was higher and 1.89 the odds of acute respiratory failure compared to the person receiving LVRS. The confidence interval here crosses 1

18.1 Pneumothorax ptx

Code
robust_ptx <- svyglm(ptx ~ vol_reduction + linps, design=robust.design, family=quasibinomial())

And then the Odds Ratio with the 95% confidence interval

Code
kable(tidy(robust_ptx, conf.int = TRUE, exponentiate = TRUE, conf.level = 0.95), digits = 2) |> kable_material()
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 0.83 0.23 -0.83 0.41 0.53 1.29
vol_reduction 0.38 0.34 -2.79 0.01 0.20 0.75
linps 0.86 0.12 -1.30 0.20 0.68 1.08

The doubly-robust odds ratio estimate for the `BLVR group was 0.38, with a 95% confidence interval (0.2,0.75). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having pneumothorax for the person receiving BLVR was lower and was 0.38 the odds of pneumothorax of the person receiving LVRS. The confidence interval did not cross 1

19 Final comparison of Outcome Estimates

19.1 Acute Respiratory Failure arf

Code
arf_blvr <- tibble(
    analysis = c("Unadjusted", "Matched", 
                 "ATT Weighted", "Double-Robust"),
    estimate = c(1.65, 1.73, 1.68, 1.89),
    conf.low = c(1, 0.82, 0.83, 0.93),
    conf.high = c(2.74, 3.63, 3.4, 3.87))

ggplot(arf_blvr, aes(x = analysis, y = estimate)) +
    geom_errorbar(aes(ymax = conf.high, ymin = conf.low), width = 0.5, color = "orange3") + 
    geom_label(aes(label = estimate), size = 6, color = "white", fill = "lightblue3") +
  labs(title = "Comparison Of Odds Ratio Estimates", subtitle = "Between Different Adjustment Approaches In Respiratory Failure", x = "Method of PS Adjustment", y = "Odds Ratio")

Method Estimate Lower 95% CI upper 95% CI
Unadjusted 1.65 1.00 2.74
Matched 1.73 0.82 3.63
Weighted 1.68 0.83 3.4
Double-Robust 1.89 0.93 3.87

Based on the comparison above, the odds ratio of acute respiratory failure didn’t substantially change from its estimate prior to propensity score adjustment. All of the confidence intervals contain 1.

19.2 Pneumothorax ptx

Code
ptx_blvr <- tibble(
    analysis = c("Unadjusted", "Matched", 
                 "ATT Weighted", "Double-Robust"),
    estimate = c(0.4, 0.68, 0.36, 0.38),
    conf.low = c(0.27, 0.38, 0.19, 0.2),
    conf.high = c(0.58, 1.22, 0.67, 0.75))

ggplot(ptx_blvr, aes(x = analysis, y = estimate)) +
    geom_errorbar(aes(ymax = conf.high, ymin = conf.low), width = 0.5, color = "purple3") + 
    geom_label(aes(label = estimate), size = 6, color = "white", fill = "lightpink3") +
  labs(title = "Comparison Of Odds Ratio Estimates", subtitle = "Between Different Adjustment Approaches In Pneumothorax", x = "Method of PS Adjustment", y = "Odds Ratio")

Method Estimate Lower 95% CI upper 95% CI
Unadjusted 0.4 0.27 0.58
Matched 0.68 0.38 1.22
Weighted 0.36 0.19 0.67
Double-Robust 0.38 0.2 0.75

The comparison above showed that the odds ratio for pneumothorax also didn’t substantially change after propensity score adjustment, reflecting a plausible conclusion that patients undergoing lung volume reduction surgery tended to have higher odds of pneumothorax based on the nature of invasive surgical intervention on the lung.

20 Consclusions

20.1 Answering The Research Questions

Is bronchoscopic lung volume reduction associated with higher odds of developing post-procedural acute respiratory failure in adult patients with emphysema compared to lung volume reduction surgery?

Yes. Based on the odds ratio obtained prior to propensity score adjustment and after, the patients undergoing bronchoscopic lung volume reduction surgery tended to have higher odds of developing acute respiratory failure, granted that the confidence interval for this estimate always crossed 1

It was noted that despite propensity score matching, weighting and a double robust adjustment (Section 19.1), the odds ratio of acute respiratory failure did not show a meaningful change from prior to adjustment. The small fluctuations wouldn’t have significant clinical implication.

Even though the confidence intervals crossed 1, this analysis showed an interesting trend in the odds of respiratory failure. A possible explanation could be that there was a higher percentage of patients with chronic respiratory failure undergoing bronchoscopic lung volume reduction in comparison to patients undergoing surgical lung volume reduction, so the patient’s lungs have already been compromised and further interventions

Is bronchoscopic lung volume reduction associated with a higher odds of developing post-procedural pneumothorax in adult patients with emphysema compared to lung volume reduction surgery?

No. BLVR was associated with lower odds of pneumothorax (Section 19.2) despite propensity score adjustment methods, which was rather plausible, as the surgical procedure was more invasive and expected to cause a higher incidence of pneumothorax. The presence of pneumothorax is well described in the literature, and could indicate a favorable response after lung volume reduction

20.2 Future Directions

Code
ggplot(blvr, aes(x = year, fill = vol_reduction_f)) +
    geom_bar(position = "dodge", col = "white") +
    scale_fill_viridis_d(option = "turbo", alpha = 0.5) +
  labs(title = "Lung Volume Reduction", subtitle = "Distribution By Year", x = "Volume Reduction", y = "# of Patients", fill = "Volume Reduction") 

Bronchoscopic lung volume reduction with endobronchial valves has largely replaced lung volume reduction surgery as a possible solution to improve excessive dyspnea. Patient selection is a key element in developing post procedural complications that could be life threatening.

Understanding the potential complications that may arise in real-world settings outside controlled clinical trials becomes increasingly crucial as the number of procedures grows and their availability expands beyond large referral centers. This knowledge aids physicians in better patient selection for the procedure, enabling them to implement effective action plans. For instance, the introduction of a “Pneumothorax Travel Cart” containing procedural kits for pneumothorax management has significantly improved time to management of pneumothorax after bronchoscopic lung volume reduction.

By trying to get a sense of the complication rates nationally, this information could be used to obtain further data from endobronchial valve device companies. These data could include more patient covariates that play a role in procedural selection, and provide a better characterization of the complications encountered.

20.3 Learning Curve

By doing this project, I got a chance to apply all the matching methods (Section 14) we learned over the span of this course. While I was slightly distraught at the beginning when I ended up with a sample with more treated patients than controls, I used it as an opportunity to go over all the plausible methods to properly troubleshoot and be prepared for similar scenarios in the future beyond the security of this class.

21 Acknowledgments

  • Dr. Marc Georgi, a gastroenterologist at Inova Fairfax Hospital, for simplifying the use of the Healthcare Utilization Project databases (the national inpatient sample and the national readmission database).
  • Dr. Thomas Love, for making statistics incredibly enjoyable!

22 References

  1. Global strategy for the diagnosis, management and prevention of chronic obstructive pulmonary disease. Global Initiative for Chronic Obstructive Lung Disease; 2021
  2. Fishman A, Martinez F, Naunheim K, et al. A randomized trial comparing lung-volume-reduction surgery with medical therapy for severe emphysema. N Engl J Med. 2003;348(21):2059-2073
  3. Fernandez-Bussy S, Labarca G, Herth FJF. Bronchoscopic Lung Volume Reduction in Patients with Severe Emphysema. Semin Respir Crit Care Med. 2018;39(6):685-692.
  4. Criner GJ, Sue R, Wright S, et al. A Multicenter Randomized Controlled Trial of Zephyr Endobronchial Valve Treatment in Heterogeneous Emphysema (LIBERATE). Am J Respir Crit Care Med. 2018;198(9):1151-1164.
  5. Criner GJ, Delage A, Voelker K, et al. Improving Lung Function in Severe Heterogenous Emphysema with the Spiration Valve System (EMPROVE). A Multicenter, Open-Label Randomized Controlled Clinical Trial. Am J Respir Crit Care Med. 2019;200(11):1354-1362
  6. Zantah M, Gangemi AJ, Criner GJ. Bronchoscopic lung volume reduction: status quo. Annals of translational medicine 2020; 8:1469.
  7. Davey C, Zoumot Z, Jordan S, McNulty WH, Carr DH, Hind MD, Hansell DM, Rubens MB, Banya W, Polkey MI, et al. Bronchoscopic lung volume reduction with endobronchial valves for patients with heterogeneous emphysema and intact interlobar fissures (the BeLieVeR-HIFi study): a randomised controlled trial. Lancet 2015;386:1066–1073
  8. Valipour A, Slebos DJ, Herth F, Darwiche K, Wagner M, Ficker JH, Petermann C, Hubner RH, Stanzel F, Eberhardt R; IMPACT Study Team. Endobronchial valve therapy in patients with homogeneous emphysema: results from the IMPACT study. Am J Respir Crit Care Med 2016;194:1073–1082.
  9. Fiorelli A, D’Andrilli A, Bezzi M, et al. Complications related to endoscopic lung volume reduction for emphysema with endobronchial valves: results of a multicenter study. J Thorac Dis. 2018;10(Suppl 27):S3315-S3325. doi:10.21037/jtd.2018.06.69
  10. Posthuma R, Vaes AW, Walraven KHM, et al. Implementation of Bronchoscopic Lung Volume Reduction Using One-Way Endobronchial Valves: A Retrospective Single-Centre Cohort Study. Respiration. 2022;101(5):476-484. doi:10.1159/000520885
  11. Healthcare Cost and Utilization Project. HCUP National Inpatient Sample (NIS). Agency for Healthcare Research and Quality; 2019-2020.

23 Session Information

Code
xfun::session_info()
R version 4.3.2 (2023-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22631)

Locale:
  LC_COLLATE=English_United States.utf8 
  LC_CTYPE=English_United States.utf8   
  LC_MONETARY=English_United States.utf8
  LC_NUMERIC=C                          
  LC_TIME=English_United States.utf8    

Package version:
  abind_1.4-5         askpass_1.2.0       backports_1.4.1    
  base64enc_0.1-3     bit_4.0.5           bit64_4.0.5        
  blob_1.2.4          boot_1.3-28.1       brio_1.1.4         
  broom_1.0.5         bslib_0.6.1         cachem_1.0.8       
  callr_3.7.3         car_3.1-2           carData_3.0-5      
  cellranger_1.1.0    chk_0.9.1           class_7.3-22       
  cli_3.6.2           clipr_0.8.0         cobalt_4.5.3       
  codetools_0.2-19    colorspace_2.1-0    compiler_4.3.2     
  conflicted_1.2.0    cpp11_0.4.7         crayon_1.5.2       
  curl_5.2.0          data.table_1.14.10  DBI_1.2.1          
  dbplyr_2.4.0        deldir_2.0-2        desc_1.4.3         
  diffobj_0.3.5       digest_0.6.34       dplyr_1.1.4        
  dtplyr_1.3.1        e1071_1.7-14        ellipsis_0.3.2     
  evaluate_0.23       fansi_1.0.6         farver_2.1.1       
  fastmap_1.1.1       fontawesome_0.5.2   forcats_1.0.0      
  foreach_1.5.2       fs_1.6.3            gargle_1.5.2       
  gbm_2.1.9           gdata_3.0.0         generics_0.1.3     
  ggformula_0.12.0    ggplot2_3.4.4       ggridges_0.5.5     
  glmnet_4.1-8        glue_1.7.0          gmodels_2.18.1.1   
  googledrive_2.1.1   googlesheets4_1.1.1 graphics_4.3.2     
  grDevices_4.3.2     grid_4.3.2          gridExtra_2.3      
  gtable_0.3.4        gtools_3.9.5        haven_2.5.4        
  highr_0.10          hms_1.1.3           htmltools_0.5.7    
  htmlwidgets_1.6.4   httr_1.4.7          ids_1.0.1          
  interp_1.1-6        isoband_0.2.7       iterators_1.0.14   
  janitor_2.2.0       jomo_2.7-6          jpeg_0.1-10        
  jquerylib_0.1.4     jsonlite_1.8.8      kableExtra_1.3.4   
  knitr_1.45          labeling_0.4.3      labelled_2.12.0    
  lattice_0.22-5      latticeExtra_0.6-30 lifecycle_1.0.4    
  lme4_1.1-35.1       lubridate_1.9.3     magrittr_2.0.3     
  MASS_7.3-60.0.1     Matching_4.10-14    Matrix_1.6-1.1     
  MatrixModels_0.5-3  memoise_2.0.1       methods_4.3.2      
  mgcv_1.9.0          mice_3.16.0         mime_0.12          
  minqa_1.2.6         mitml_0.4-5         mitools_2.4        
  modelr_0.1.11       mosaic_1.9.0        mosaicCore_0.9.4.0 
  mosaicData_0.20.4   munsell_0.5.0       naniar_1.0.0       
  nlme_3.1-163        nloptr_2.0.3        nnet_7.3-19        
  norm_1.0.11.1       numDeriv_2016.8.1.1 openssl_2.1.1      
  ordinal_2023.12.4   pan_1.9             parallel_4.3.2     
  pbkrtest_0.5.2      pillar_1.9.0        pkgbuild_1.4.3     
  pkgconfig_2.0.3     pkgload_1.3.3       plyr_1.8.9         
  png_0.1-8           praise_1.0.0        prettyunits_1.2.0  
  processx_3.8.3      progress_1.2.3      proxy_0.4-27       
  ps_1.7.5            purrr_1.0.2         quantreg_5.97      
  R6_2.5.1            ragg_1.2.7          rappdirs_0.3.3     
  RColorBrewer_1.1-3  Rcpp_1.0.12         RcppEigen_0.3.4.0.0
  readr_2.1.5         readxl_1.4.3        rematch_2.0.0      
  rematch2_2.1.2      repr_1.1.6          reprex_2.1.0       
  rlang_1.1.3         rmarkdown_2.25      rpart_4.1.23       
  rprojroot_2.0.4     rstudioapi_0.15.0   rvest_1.0.3        
  sass_0.4.8          scales_1.3.0        selectr_0.4.2      
  shape_1.4.6         skimr_2.1.5         snakecase_0.11.1   
  SparseM_1.81        splines_4.3.2       stats_4.3.2        
  stringi_1.8.3       stringr_1.5.1       survey_4.2-1       
  survival_3.5-7      svglite_2.1.3       sys_3.4.2          
  systemfonts_1.0.5   tableone_0.13.2     testthat_3.2.1     
  textshaping_0.3.7   tibble_3.2.1        tidyr_1.3.0        
  tidyselect_1.2.0    tidyverse_2.0.0     timechange_0.2.0   
  tinytex_0.49        tools_4.3.2         twang_2.6          
  tzdb_0.4.0          ucminf_1.2.1        UpSetR_1.4.0       
  utf8_1.2.4          utils_4.3.2         uuid_1.1.1         
  vctrs_0.6.5         viridis_0.6.4       viridisLite_0.4.2  
  visdat_0.6.0        vroom_1.6.5         waldo_0.5.2        
  webshot_0.5.5       withr_2.5.2         xfun_0.41          
  xgboost_1.7.7.1     xml2_1.3.6          xtable_1.8-4       
  yaml_2.3.8          zoo_1.8-12